Diffusion in a crowded environment 
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We analyze a pair of diffusion equations which are derived in the infinite system-size limit from 
a microscopic, individual-based, stochastic model. Deviations from the conventional Fickian pic- 
ture are found which ultimately relate to the depletion of resources on which the particles rely. 
The macroscopic equations are studied both analytically and numerically, and are shown to yield 
anomalous diffusion which does not follow a power law with time, as is frequently assumed when 
fitting data for such phenomena. These anomalies are here understood within a consistent dynam- 
ical picture which applies to a wide range of physical and biological systems, underlining the need 
for clearly defined mechanisms which are systematically analyzed to give definite predictions. 
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Almost all discussions relating to the modeling of diffu- 
sion assume Pick's law — that the rate at which one sub- 
stance diffuses through another is directly proportional to 
the concentration gradient of the diffusing substance [l| . 
There is a good reason for this: it is empirically very well 
supported, at least at not too high concentrations, and is 
part of the wider theoretical framework of linear nonequi- 
librium thermodynamics However one might expect 
nonlinear corrections in various situations e.g. if obstacles 
arc present or if there is more than one species at high 
concentrations, and indeed many reports of anomalous 
behavior can be found in the literature These 
studies are mainly experimental, or involve numerical 
simulations, and the theory which is given mainly con- 
sists of phenomenological fits to the data. Many of the 
fits suggest that the mean-square displacement of the dif- 
fusing species grows with time, t, like t", where a < 1. 
The phenomena described go under various names such 
as "molecular crowding" |12h14| and "single-file diffu- 
sion" [TB - flsl , all indicating that the increased concen- 
tration impedes the flow of particles in some way, leading 
to a violation of Pick's law. 

There are relatively few first-principles studies of these 
effects. Those that do exist include hydrodynamical 
models with nonlinear constitutive relations [Toj and 
simple symmetric exclusion processes in one dimen- 
sion [l^. In this Letter we propose a modification of 
Pick's law which is based on a physically motivated mi- 
croscopic theory. Unlike previous theories of single-file 
diffusion [2^ [2l| it holds in arbitrary dimensions, since 
it is not a consequence of the physical 'jamming' of the 
particles, rather it is due to the depiction of resources 
on which the particles rely. This resource may be space 
to move, but could also be a chemical substrate required 
for the system to remain viable. The phenomenon we 
describe is seen in experiments which trace the motion 
of tagged particles. This is equivalent to assuming the 
existence of two types of particle with the same diffusion 
constant. 

We begin with a generic microscopic system in a given 
volume of d— dimensional space which is divided into a 



large number, £7, of small (hypercubic) patches. Pach 
patch, labeled by a, can contain up to N particles: tIq of 
type A, rria of type B, and = N — — rua vacancies, 
denoted by E. We will assume that the particles have no 
direct interaction, however there will be an indirect in- 
teraction in that the mobility of particles will be affected 
if neighboring patches have few vacancies. 

To model this more concretely we assume that the par- 
ticles move only to nearest-neighbor patches, and then 
only if there is a vacancy there: 



A 



/3: 



Ba + Ep Ea + Bp. 



(1) 



Here a and /3 label nearest- neighbor patches with Aa , Ba , 
and Ea being the respective types of particles in patch 
a, and /ii and /i2 being the reaction rates. The state 
of the system will be characterized by the number of 
A and B particles in each patch, that is, by the vec- 
tor n = (rii, . . . , nji), where ria = (tIq, ma). The rate of 
transition from state n', to another state n, is denoted by 
T(n|n') — with the initial state being on the right. The 
transition rates associated with the migration between 
nearest-neighbor patches take the form 



T{na - I, Tip + l\na,np) 
T{ma — 1, mp + l\ma,mp) 



/ii Ua N — np — rap 

'Tnlv N ' 

H2 mg N - np - mp 



(2) 



where z is the number of nearest neighbors that each 
patch has and where within the brackets we have chosen 
to indicate only the dependence on those particles which 
are involved in the reaction. It is the presence of the 
factor vp — N ~ np ~ mp, which reduces the transition 
rate if there are few vacancies in the target patch, and 
which modifies Pick's law in the macroscopic theory. 

This is a Markov process, and the probability of finding 
the system in state n at time t, denoted by P(n, t), is 



2 



given by the master equation 



dP(n,t) 
di 



J2 [T{n\n')P{n',t)-T{n'\n)P{n,t)] 



(3) 

where the aUowed transitions are those given by Eq. ([T]|. 
This defines the microscopic process, but we are inter- 
ested in the macroscopic equations that this process gen- 
erates. To find these we need to find the dynamical equa- 
tions for the ensemble averages (ua) and {nia). Multi- 
plying Eq. ([3]) by n„ and summing over all n gives, after 
shifting some of the sums by ±1, 



dt 



[{T{na + 1, n/3 - l\na,ni3)) 



pea 



{T{na -l,np + l|nQ,n^))] , 



(4) 



where the notation X),9eQ ni^ans 'sum over all patches /3 
which are nearest-neighbors of the patch a\ A similar 
equation holds for d{ma)/dt. 

The averages in Eq. (|4]) are carried out by using the 
explicit forms ^ and replacing the averages of prod- 
ucts by the products of averages, which is valid in the 
limit N ^ oo. Then scaling time by a factor of Nfl one 
finds ^22] 



dt 
dtpa 
dt 



Here 



ba = hm ipa = lim 



(6) 



and A is the discrete Laplacian operator defined by 
A/a ^ (2/2:)E/36q(//3 - fa)- Finally taking the size 
of the patches to zero, and scaling the rates fii and /i2 
appropriately [22| to give diffusion constants Di and D2 , 
gives partial differential equations for (/)(x, t) and t/'(x, t): 

^ = D2 [V2^ + i'V^<P ~ 0VV] , (7) 

where is the usual Laplacian. 

We can give a quite complete analysis of Eqs. in the 
case when the diffusion constants are equal. Let Di — 
D2 = D and absorb D into the definition of the time. 
Then adding the two equations gives 



(8) 



whereas the equation for the difference a = {(f) ~ 4')/V^ 
is 

da 



dt 



a/2 (ctV V - pV V) 



(9) 



We will take initial conditions such that p(— x, 0) = 
p(x,0) and cr(-x,0) = -cr(x,0). Solving Eq. dH) for 
p and going over to Fourier space gives for Eq. (|9|): 



(10) 



da{k,t) 

dt ' ■ (27r) 

xp(k-p,0)e-('^-P)'*a(p,t), 



where p(k — p, 0) is the initial value of p in Fourier space. 
This equation is linear in a which allows us to make fur- 
ther analytic progress. We have considered two types of 
particle for simplicity; the above discussion can easily be 
extended to three or more types. 

We will first analyze Eq. ([TU]) in one dimension. The 
calculation in d— dimensions is not much more difficult, 
and we will give the main result later, but is less clear 
due to the number of indices involved in the interme- 
diate steps. We begin by noting that p{k,t) is even 
in k and a{k,t) is odd in k, so that cr(0,i) = 0. So 
the first non-trivial term in the expansion of a{k,t) is 
{x{t))a = —id(j{k,t)dk\k=o- From Eq. (fTO|) one finds 
that 



Ax{t))c 

' dt 



^2V2j ^pp{p,0)e-p''a{p,t). (11) 



In order to evaluate the integral we need to know some- 
thing of the behavior of a{p,t). Numerical simula- 
tions (sec later) indicate that the ratio of (a;^"+^(t))o- to 
{t))a, n = 1, 2, . . ., is proportional to i, for large t. 



J2n-1 



(5) to a very good approximation. Therefore 



a{p,t) 



n=0 



2" + l(x2"+l(t)), 



(2n + l)! 



i{x{t))„pf{pH), 



for large t where 



n=0 



(-i)"a«iy" 

(2n-t-l)! ' 



(12) 



(13) 



and where the a„ are constants. Substituting Eq. (jl2p 
into Eq. (fTTj) and scaling p^ by t gives 



d{x{t))a Ai 



where 



dt 2i3/2 



AVi 



^i = — / dpf{p')p'e 



(14) 



(15) 



This scaling implies that p{p, 0) is replaced by p(j) = 
0,0) = / dxp{x,0) = V2V1, for large t. Solving the dif- 
ferential equation (fT4| gives 



{x{t))a = cexp 



VtJ' 



(16) 



where c is a constant. 
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The coefficients a„ in Eq. p3)) may be found by differ- 
entiating Eq. pU)) {2n + 1) times with respect to k and 
then setting fc = 0. In the resulting differential equation 
for {x^"-~^^ (t)) use of Eq. shows that the p— integral 
is down by powers of t on the contribution coming from 
the term —k'^a{k,t) for n > 0. So for large t, 



dt 



2n(2n+l)(x2"-i(t)), 
2n(2n + l)a„ 



(17) 



Using Eq. p6)) the integration in Eq. ([TT]) can be car- 



ried out for large t. This gives 



2(2n -t- 

l)a„_it"(a;(t))o. and so a„ = 2(2ri, + l)a„_i. There- 
fore from Eq. ([T^ /(y) and from Eq. Ai ~ 
Vi/V27r. Finally, from Eqs. ^ and (HH) we find that 

cr(fc,t) ^ ifccexp ^— — ^=^^ cxp (— fc^t), (18) 

for large t. Taking the inverse Fourier transform gives 



a{x, t) 



cx 



47^1/2^3/2 



exp 



cxp(-2;V4t). (19) 



The corresponding calculation in d— dimensions can in 
principle be carried out in a similar way. This will be 
discussed in more detail elsewhere and here we only give 
the generalization of Eq. ()16p : 



Ci exp 



'i-d/2 



(20) 



where the Ci, 
by 



1, . . . , d, are constants and is given 



dp 



Vd 



2{3d~2)/2^d/2^- 



(21) 



Returning to the one-dimensional case, we may use the 
results we have obtained to characterize the time evolu- 
tion of the variance of the original distributions (f) and 
ijj. Since p has only even moments and cr only odd mo- 
ments, {x^'')cii = (— l)"(a;")^, for any integer n. It follows 
that {{x^))^ = {{x^))^, where the symbol ((•))/ stands 
for the variance of function /. Using the result for 
the first moment and the standard diffusion result for the 
second moment, we have for the normalized mean-square 
displacement: 



((a;'(t)))0> = 2<. 



2V^ 



exp 



1/2 



V2Vi 

(22) 

For large enough times, the system displays normal dif- 
fusion: the variances scale linearly with time, and are 
shifted by a constant factor. For relatively short times, 
of the order of the inverse of the diffusion coefficient, here 
absorbed in the definition of t, deviations from the usual 
behavior are predicted to occur due to the exponential 
factor in Eq. (P^ . which reduces the diffusion. 




FIG. 1: (Color online) Main panel: the variances {{x))^ and 
{{x))^ are plotted as a function of the rescaled time, Dt. Sym- 
bols refer to the numerical simulations, diamonds and pluses 
stand respectively for (j> and ip. The dashed line represents the 
normal diffusion prediction, while the thick solid line refers to 
formula ^ where c = 2.035. Here D = 0.003 and Vi = 1. 
Inset: The first moment {x)a is represented as a function of 
rescaled time, Dt. Symbols refer respectively to D = 0.003 
(circles) and D = 0.008 (triangles). The solid line is the 
asymptotic prediction with c = 2.035. 



The validity of the approximations we have made in 
the above analysis have been checked by carrying out 
numerical simulations for the one-dimensional case. We 
used Euler discretization, both in the space and time co- 
ordinates of Eq. ([7|) and assumed identical diffusion con- 
stants for both species, so as to make contact with the 
analysis. In the main panel of Fig. [T] the time evolution 
of the measured variances, for both <j) and tp, arc dis- 
played. The variances are seen to grow more slowly than 




FIG. 2: (Color online) Left panel: three profiles of the distri- 
bution (f> are plotted (thin line), corresponding to successive 
instants of the dynamics. From top to bottom, the snapshots 
are taken at rescaled time Dt = 0.3,0.9, 1.5. The theoreti- 
cal profile at Dt — 1.5 is represented with a thick line (red 
online). Right panel: Same as in left panel, but for the dis- 
tribution if). The analytical prediction is plotted with a thick 
line (blue online). The simulations refer to D = 0.003 and 
Vi = 1. 
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FIG. 3: (Color online) The variances {{x))^ (diamonds) and 
((x))^ (pluses) are plotted as a function of their respective 
rescaled time, Dat, with a = 1, 2. In both cases the variances 
grow slower than predicted by the standard diffusion theory 
(dashed line). Here, Di = 0.003, D2 = 0.008 and Vi = 1. 



predicted from standard diffusion theory. The agreement 
between the theoretical prediction ([^^ and the simula- 
tions is excellent, even at quite short times. The constant 
c is determined by fitting the late time evolution of {x)a 
from simulations to the asymptotic profile (|16p . see in- 
set of Fig. [TJ The form of the time evolution of the first 
moment {x)a given by Eq. (fT6)) has also been checked 
numerically for a range of parameters, and is always in 
excellent agreement. As a further check. Fig. [2] shows 4> 
and -0 found from numerical simulations for a range of 
times. The recorded snapshots at the largest time show 
good agreement with the theoretical curves at that time. 

These considerations all point to the validity of the 
analysis for the case when the diffusivities of the two 
species are equal. To shed light on the more general set- 
ting where D\ ^ I?2 7 and so broaden the range of applica- 
bility of our conclusions, we rely on numerical simulation. 
In Fig. [3l we show the time evolution of the variance of 



(resp. ■(/;) vs. the rescaled time D\t (resp. D^l). As 
in the symmetric case, the growth of the variances is 
slower than for standard diffusion. This again reflects 
the presence of the finite carrying capacity imposed at 
the level of the microscopic dynamics. Remarkably, the 
function {x)a is again found to empirically obey Eq. (|16p . 
the undetermined factor A\ depending on the individual 
diffusivities. 

The modified diffusive behavior we have found is de- 
rived from a general principle formulated at the micro- 
scopic level, not from a plienomenological fit. It is im- 
portant to stress this fact, since virtually all the work on 
molecular crowding, and related phenomena, to date has 
postulated that the mean-square displacement increased 
in time like a power: ((x^)) ~ i". There has been much 
discussion of whether a < 1 (called 'subdiffusion') or 
a > 1 (called 'superdiffusion'). Our results can be fit- 
ted by either, depending on when the time-window for 
the fit is taken; if taken at early times subdiffusion is 
found, at late times superdiffusion is found. Indeed by 
picking suitable time-windows a wide range of values of 
a can be found. This ambiguity shows the necessity of 
starting with a clearly defined mechanism, which can be 
precisely implemented (not phcnomenologically invoked) 
and from which clear predictions can be systematically 
derived. The approach we have described, and the re- 
sults we have obtained, in this Letter follow this philoso- 
phy closely and we expect that comparison of our results 
with experiments in the future will help to clarify the ef- 
fect of molecular crowding and of resource depletion on 
diffusion. 
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